import os
from array import *
import csv
import math
import random
import numpy as np
import scipy as sp
import scipy.stats.stats as st
from scipy.interpolate import UnivariateSpline
from scipy.optimize import curve_fit
import scipy.stats.stats as st
from scipy.stats import t,norm,triang
#from lmfit import minimize, Parameters,report_errors

import matplotlib.pyplot as plt
import matplotlib.mlab as mlab

from matplotlib.ticker import MaxNLocator, AutoMinorLocator
from matplotlib import transforms as mtransforms
from matplotlib import font_manager as font
from matplotlib.figure import Figure #as mFig
from matplotlib.axis import Axis




#constants
E=0.735        #Nm
x1=0.073       #m
k=2*E/(x1**2)
x=0.005        #m
g=9.81         #m/s**2
M=k*x/g     #kg
V1sq=2*E/M   #(m/s)**2


def convert_Rt_to_Rh(Rt,theta):
    x2t=Rt*x1/100.
    V2tsq=(1/(0.5*M))*(0.5*k*x2t**2 + M*g*x2t*np.cos(np.radians(theta)))
    V2sq=V2tsq*0.5*k*x1**2/(0.5*k*x1**2+M*g*x1*np.cos(np.radians(theta)))
    Rh=np.sqrt(V2sq/V1sq)*100
    return Rh

theta=0,30,60,90,120
Rt=30

print theta
print convert_Rt_to_Rh(Rt,theta)

def SMD_analysis(x,numsamples,samp_size,alpha):
    #creating sample mean distribution
    cur_smd=samplemeandist(x,numsamples,samp_size)

    #checking if sample means are normally distributed
    pKS=KS2samp_normtest(cur_smd,alpha)
    
    #6.3 plotting fitted normal distribution curve
    mu, std = norm.fit(cur_smd)
    print "mean: ",round(mu,0), "   stdev: ", round(std,0)
    print round(mu-2*std,0), round(mu+2*std,0)
    #plt.hist(cur_smd, bins=25, normed=True, alpha=0.6,histtype='step',color='0.6')
    x = np.linspace(min(cur_smd),max(cur_smd), 200)
    p = norm.pdf(x, mu, std)
    plt.plot(x, p, '-',linewidth=2, label=lith[f])
    
    return cur_smd

def samplemeandist(x,numsamples,samp_size):
    sample_means=array('f')
    for a in range(numsamples):
        indices = np.arange(len(x))
        np.random.shuffle(indices)
        samp_indices=indices[:samp_size]
        sample=x[samp_indices]
        sample_means.append(np.mean(sample))
    return sample_means

def KS2samp_normtest(x,alpha):
    normal_samp=np.random.normal(np.mean(x),np.std(x),len(x))
    KS,pKS=st.ks_2samp(x, normal_samp)
    if pKS>alpha:print "dist is normal; KStest pval=", round(pKS,2)
    else:print "dist is NOT normal; KStest pval=", round(pKS,2)
    return pKS


np.random.seed(4657)


try:
    fo1 = open("sleyteschmidthammerdata2014.csv", 'rb')
    print "file found"
    fo = csv.reader(fo1,delimiter=',')
    lines_appended=0
    lithgroup=[]
    R=[]
    theta=[]
    rcount=0
    for i in fo:
        try:
            if float(i[9])<10.:
                print i[6],i[9]
                continue
            rcount=rcount+1
            R.append(float(i[9]))
            #print min(R),max(R)
            theta.append(float(i[8]))
            if i[6]=="Cce" or i[6]=="Ccw":
                lithgroup.append("CC")
            else:
                lithgroup.append(i[6])
            #print i[6],i[8],i[9]
        except:
            continue
    print "Rread: ",rcount
except:
    print "missing file"


#generating function for Rh-to-UCS conversion
fo1 = open("RtoUCS.dat", 'rb')
fo = csv.reader(fo1,delimiter=' ')
R2=[]
UCS=[]
R2.append(10)
R2.append(0)
UCS.append(5)
UCS.append(1)
for i in fo:
    R2.append(float(i[0]))
    UCS.append(float(i[2]))
R2=np.asarray(R2)
UCS=np.asarray(UCS)
par=np.polyfit(R2,UCS,3)
func=np.poly1d(par)



print lithgroup


go=1
if go==1:
    lith=np.asarray(['CC','FC','LS','UM'])

    smd=[]
    fig,axes=plt.subplots(nrows=len(lith),ncols=1,sharex=True,sharey=True,figsize=(5,7))
    fig2,axes2=plt.subplots(nrows=len(lith),ncols=1,sharex=True,sharey=True, figsize=(4,7))
    fig3,axes3=plt.subplots(nrows=len(lith),ncols=1,sharex=True,sharey=True,figsize=(4,7))
    for f in range(len(lith)):
        
        i=np.where((np.asarray(lithgroup)==lith[f]))
        curR=np.asarray([R[a] for a in i[0]])                      
        curtheta=np.asarray([theta[a] for a in i[0]])
        curRh=convert_Rt_to_Rh(curR,curtheta)
        print len(curRh)

        curax=axes[f]
        plt.sca(curax)
        plt.title("\n"+lith[f]+" ,n="+str(len(curRh)),fontsize='small')
        curax.hist(curRh,bins=int((max(curRh)-min(curRh))/2),normed=True)
        if f==len(lith)-1:
            curax.set_xlabel("Rebound number (horizontal)")
        plt.setp(curax.get_yticklabels(), visible=False)
        plt.tight_layout()


        
        curUCS=func(curRh)
        
        curax=axes2[f]
        plt.sca(curax)
        plt.title("\n"+lith[f]+" ,n="+str(len(curRh)),fontsize='small')
        
        curax.hist(curUCS,bins=int((max(curUCS)-min(curUCS))/2),normed=True)

        print "*******************************"
        print lith[f]
        params=triang.fit(curUCS)
        print params[1], params[1]+params[0]*params[2], params[1]+params[2]
        h=1*2/params[2]
        curax.plot([params[1], params[1]+params[0]*params[2], params[1]+params[2]],[0,h,0],'r-')

        KS2samp_normtest(curUCS,0.05)

        if f==len(lith)-1:
            curax.set_xlabel("UCS, MPa")
        plt.setp(curax.get_yticklabels(), visible=False)
        plt.tight_layout()

        
       
        
        curax=axes3[f]
        plt.sca(curax)
        plt.title("\n"+lith[f],fontsize='small')
        numsamples=500
        samp_size=10
        alpha=0.05
        print "\n",lith[f],round(np.mean(curRh),0)
        smd.append(SMD_analysis(curUCS,numsamples,samp_size,alpha))
        curax=plt.gca()
        if f==len(lith)-1:
            curax.set_xlabel("UCS, MPa")
        plt.setp(curax.get_yticklabels(), visible=False)
        plt.tight_layout()

    

    for f in range(len(lith)):
        cursmd=smd[f]
        print ""
        for g in range(f+1,len(lith)):
            comsmd=smd[g]

            T,pT=st.ttest_ind(cursmd,comsmd)
            KS,pKS=st.ks_2samp(cursmd,comsmd)
            print [lith[f], lith[g], round(pT,3), round(pKS,3)]
                          
                          
plt.show()        
        
                
                
            
            
